library(rms)
library(tidyverse)
library(corrplot)
library(dplyr)
library(RColorBrewer)
library(plotly)
library(sjPlot)
require(foreign)
require(MASS)
require(Hmisc)
require(reshape2)
require(plm)
library(lmtest)
library(ggthemes)
library(tidyr)
require(ggeffects)
library(ordinal)
library(broom)
require(lmtest) 
library(stargazer)
library("sandwich")
require(corrgram)
library(data.table)



# R Version 2023.06.2+561 (2023.06.2+561)
# setwd
dat <- read.csv("spotlight_final.csv")


########################################
#################### PR*peer review
########################################

dat$outcome_final <- as.factor(dat$outcome_final)

m_1 <- polr(outcome_final ~  
              peerrev + pr + startyear, data = dat, Hess = TRUE)
summary(m_1)
m_1 <- coeftest(m_1, vcov=vcovCL(m_1, factor(dat$NCP)))
m_1


m_2 <- polr(outcome_final ~  
              peerrev + pr + NGOnumbers + incumpower + opppower +
              provision + provision2 + startyear, data = dat, Hess = TRUE)


ggpredict(m_2, term = c("peerrev"), 
          condition = c(startyear = mean(dat$startyear),
                        pr = 1, 
                        provision = mean(dat$provision, na.rm=T), incumpower = mean(dat$incumpower, na.rm=T), 
                        opppower = mean(dat$opppower, na.rm=T),
                        NGOnumbers = mean(dat$NGOnumbers),
                        provision2 = mean(dat$provision, na.rm=T)*mean(dat$provision, na.rm=T)), 
          vcov.fun = "vcovCL", vcov.type = "HC1", 
          vcov.args = list(cluster = dat$subclass))

m_2 <- coeftest(m_2, vcov=vcovCL(m_2, factor(dat$NCP)))
m_2


m_3 <- polr(outcome_final ~  
              peerrev + pr + NGOnumbers + incumpower + opppower +
              provision + provision2 + pr*peerrev + startyear, data = dat, Hess = TRUE)
summary(m_3)
m_3 <- coeftest(m_3, vcov=vcovCL(m_3, factor(dat$NCP)))
m_3

dat$Peer_char <- as.factor(dat$Peer_char)



m_3p <- polr(outcome_final ~  
               Peer_char + PR + NGOnumbers + incumpower + opppower +
               provision + provision2 + PR*Peer_char + startyear, data = dat, Hess = TRUE)
summary(m_3p)

### Figure 3 -- peer review effects
plot(ggpredict(m_3p, term = c("Peer_char", "PR"), 
               condition = c(startyear = mean(dat$startyear),
                             provision = mean(dat$provision, na.rm=T), 
                             incumpower = mean(dat$incumpower, na.rm=T), 
                             opppower = mean(dat$opppower, na.rm=T),
                             NGOnumbers = mean(dat$NGOnumbers),
                             provision2 = mean(dat$provision, na.rm=T)*mean(dat$provision, na.rm=T)), 
               vcov.fun = "vcovCL", vcov.type = "HC1", 
               vcov.args = list(cluster = dat$subclass)), connect.lines = TRUE, ci.style = "dot") + 
  xlab("") + 
  ylab("Predicted Probabilities of \n Vindication, Neutral, Recommendation") + 
  ggtitle("Peer Review Effects") + 
  theme_classic() + scale_color_grey(start=0.7, end = 0.2)

ggpredict(m_3p, term = c("Peer_char", "PR"), 
          condition = c(startyear = mean(dat$startyear),
                        provision = mean(dat$provision, na.rm=T), 
                        incumpower = mean(dat$incumpower, na.rm=T), 
                        opppower = mean(dat$opppower, na.rm=T),
                        NGOnumbers = mean(dat$NGOnumbers),
                        provision2 = mean(dat$provision, na.rm=T)*mean(dat$provision, na.rm=T)), 
          vcov.fun = "vcovCL", vcov.type = "HC1", 
          vcov.args = list(cluster = dat$subclass))

m_3p <- coeftest(m_3p, vcov=vcovCL(m_3p, factor(dat$NCP)))
m_3p

m_4 <- polr(outcome_final ~  
              peerrev + pr  + NGOnumbers + incumpower + opppower +  treaty  + binarypta +
              provision + provision2 + pr*peerrev  + Labor + Environment + Bribery + startyear, data = dat, Hess = TRUE)

m_4 <- coeftest(m_4, vcov=vcovCL(m_4, factor(dat$NCP)))
m_4




m_4p <- polr(outcome_final ~  
               Peer_char + PR + NGOnumbers + incumpower + opppower + treaty + binarypta + 
               provision + provision2 + PR*Peer_char + Labor + Environment + Bribery + startyear, data = dat, Hess = TRUE)
summary(m_4p)

ggpredict(m_4p, term = c("Peer_char", "PR"), 
          condition = c(startyear = mean(dat$startyear),
                        provision = mean(dat$provision, na.rm=T), 
                        incumpower = mean(dat$incumpower, na.rm=T), 
                        opppower = mean(dat$opppower, na.rm=T),
                        NGOnumbers = mean(dat$NGOnumbers),
                        provision2 = mean(dat$provision, na.rm=T)*mean(dat$provision, na.rm=T)), 
          vcov.fun = "vcovCL", vcov.type = "HC1", 
          vcov.args = list(cluster = dat$subclass))

m_4p <- coeftest(m_4p, vcov=vcovCL(m_4p, factor(dat$NCP)))
m_4p


# Table 1 in the manuscript
stargazer(m_1, m_2, m_3, m_4, digits=2, style = "IO")


########################################
#################### PR*domestic NGOs
########################################
d1 <- polr(outcome_final ~  
             polipush + pr + startyear, data = dat, Hess = TRUE)

summary(d1)
d1 <- coeftest(d1, vcov=vcovCL(d1, factor(dat$NCP)))
d1



d2 <- polr(outcome_final ~  
             polipush + pr + NGOnumbers + incumpower + opppower + 
             provision + provision2 +startyear, data = dat, Hess = TRUE)
summary(d2)

ggpredict(d2, term = c("polipush"), 
          condition = c(startyear = mean(dat$startyear),
                        pr = 1,
                        provision = mean(dat$provision, na.rm=T), 
                        incumpower = mean(dat$incumpower, na.rm=T), 
                        opppower = mean(dat$opppower, na.rm=T),
                        NGOnumbers = mean(dat$NGOnumbers),
                        provision2 = mean(dat$provision, na.rm=T)*mean(dat$provision, na.rm=T)), 
          vcov.fun = "vcovCL", vcov.type = "HC1", 
          vcov.args = list(cluster = dat$subclass))

d2 <- coeftest(d2, vcov=vcovCL(d2, factor(dat$NCP)))
d2



d3 <- polr(outcome_final ~  
             polipush + pr  + NGOnumbers + incumpower + opppower +
             provision + provision2 + polipush*pr + startyear, data = dat, Hess = TRUE)
ggpredict(d3, term = c("polipush", "pr"), 
          condition = c(startyear = mean(dat$startyear),
                        provision = mean(dat$provision, na.rm=T), 
                        incumpower = mean(dat$incumpower, na.rm=T), 
                        opppower = mean(dat$opppower, na.rm=T),
                        NGOnumbers = mean(dat$NGOnumbers),
                        provision2 = mean(dat$provision, na.rm=T)*mean(dat$provision, na.rm=T)), 
          vcov.fun = "vcovCL", vcov.type = "HC1", 
          vcov.args = list(cluster = dat$subclass))


summary(d3)
d3 <- coeftest(d3, vcov=vcovCL(d3, factor(dat$NCP)))
d3


dat$homengo <- NA
for(i in 1:nrow(dat)){
  if(is.na(dat$polipush[i]) == TRUE){dat$homengo[i] <- NA}
  else if(dat$polipush[i] > 0){dat$homengo[i] <- "2.Domestic NGO"}
  else{dat$homengo[i] <- "1.Foreign"}
}

dat$Advocacy <- as.factor(dat$homengo)

d3f <- polr(outcome_final ~  
              Advocacy + PR  + NGOnumbers + incumpower + opppower +
              provision + provision2 +  startyear + PR*Advocacy, data = dat, Hess = TRUE)
summary(d3f)

## Figure 3 -- Domestic NGO Effect
plot(ggpredict(d3f, term = c("Advocacy", "PR"), 
               condition = c(startyear = mean(dat$startyear),
                             provision = mean(dat$provision, na.rm=T), 
                             incumpower = mean(dat$incumpower, na.rm=T), 
                             opppower = mean(dat$opppower, na.rm=T),
                             NGOnumbers = mean(dat$NGOnumbers),
                             provision2 = mean(dat$provision, na.rm=T)*mean(dat$provision, na.rm=T)), 
               vcov.fun = "vcovCL", vcov.type = "HC1", 
               vcov.args = list(cluster = dat$subclass)), connect.lines = TRUE, ci.style = "dot") + 
  xlab("") + 
  ylab("Predicted Probabilities of \n Vindication, Neutral, Recommendation") + 
  ggtitle("Domestic NGO Effects") + 
  theme_classic() + scale_color_grey(start=0.7, end = 0.2)



ggpredict(d3f, term = c("Advocacy", "PR"), 
          condition = c(startyear = mean(dat$startyear),
                        provision = mean(dat$provision, na.rm=T), 
                        incumpower = mean(dat$incumpower, na.rm=T), 
                        opppower = mean(dat$opppower, na.rm=T),
                        NGOnumbers = mean(dat$NGOnumbers),
                        provision2 = mean(dat$provision, na.rm=T)*mean(dat$provision, na.rm=T)), 
          vcov.fun = "vcovCL", vcov.type = "HC1", 
          vcov.args = list(cluster = dat$subclass))

d3f <- coeftest(d3f, vcov=vcovCL(d3f, factor(dat$NCP)))
d3f




d4 <- polr(outcome_final ~  
             polipush + pr + NGOnumbers + incumpower + opppower +
             provision + provision2 +  polipush*pr +  treaty + binarypta + 
             Labor + Environment + Bribery + startyear , data = dat, Hess = TRUE)

summary(d4)
d4 <- coeftest(d4, vcov=vcovCL(d4, factor(dat$NCP)))
d4


d4f <- polr(outcome_final ~  
              Advocacy + PR  + NGOnumbers + incumpower + opppower + treaty + binarypta + Labor + Environment + Bribery +
              provision + provision2 +  startyear + PR*Advocacy, data = dat, Hess = TRUE)
d4f <- coeftest(d4f, vcov=vcovCL(d4f, factor(dat$NCP)))
d4f



# homengo as a binary variable
d3_binary <- polr(outcome_final ~  
                    homengo + pr  + NGOnumbers + incumpower + opppower +
                    provision + provision2 + homengo*pr + startyear, data = dat, Hess = TRUE)
d3_binary <- coeftest(d3_binary, vcov=vcovCL(d3_binary, factor(dat$NCP)))
d3_binary


d4_binary <- polr(outcome_final ~  
                    homengo + pr + NGOnumbers + incumpower + opppower +
                    provision + provision2 +  homengo*pr +  treaty + binarypta + 
                    Labor + Environment + Bribery + startyear , data = dat, Hess = TRUE)

summary(d4_binary)
d4_binary <- coeftest(d4_binary, vcov=vcovCL(d4_binary, factor(dat$NCP)))
d4_binary

## Table 2 in the manuscript
stargazer(d1, d2, d3, d4, d3_binary, d4_binary, digits=2, style = "IO")


#######################################
############# alternative hypotheses (Appendix E. Third Party NGO)
#######################################

a11 <- polr(outcome_final ~  
              thirdpush + NGOnumbers +  startyear, data = dat, Hess = TRUE)

a11 <- coeftest(a11, vcov=vcovCL(a11, factor(dat$NCP)))
a11


a3 <- polr(outcome_final ~  
             peerrev + NGOnumbers + polipush + incumpower + 
             opppower  + pr + thirdpush + provision + provision2 + treaty +
             Labor + Environment + Bribery + startyear, data = dat, Hess = TRUE)

a3 <- coeftest(a3, vcov=vcovCL(a3, factor(dat$NCP)))
a3

## Table 7 in Appendix E
stargazer(a11, a3, digits=2, style = "IO")


#################################################
############# Matching (Appendices B and C)
#################################################

library(MatchIt)
df = na.omit(subset(dat, select = c(outcome_final, NGOnumbers, peerrev, pr, incumpower, opppower, treaty, binarypta,
                                    polipush, legparty,
                                    provision, provision2, startyear, endyear, Labor, Environment, Bribery, NCP)))


m.out0 <- matchit(peerrev ~ startyear + legparty, data = df,
                  method = "cem", cutpoints = list(startyear =2011), distance = "glm")

summary(m.out0, un = FALSE) 

# Figure 1 in Appendix B
plot(summary(m.out0), xlim=c(0,2))

############### Peer Review as a treatment in matching

mdat <- match.data(m.out0)


mdat$peerrev <- as.factor(mdat$peerrev)
mdat$pr <- as.factor(mdat$pr)

mdat$PR_char <- NA

for(i in 1:nrow(mdat)){
  if(mdat$pr[i] == 1){mdat$PR_char[i] <- "PR"}
  else if(mdat$pr[i] == 0){mdat$PR_char[i] <- "Non-PR"}
}
mdat$PR_char <- as.factor(mdat$PR_char)


for(i in 1:nrow(mdat)){
  if(mdat$peerrev[i] == 1){mdat$Peer[i] <- "Peer Review"}
  else if(mdat$peerrev[i] == 0){mdat$Peer[i] <- "No Review"}
}

mdat$Peer <- as.factor(mdat$Peer)


m_1 <- polr(outcome_final ~ Peer  +  PR_char  + startyear, data = mdat, Hess = TRUE)
m_1 <- coeftest(m_1, vcov=vcovCL(m_1, factor(mdat$NCP)))
m_1


m_2 <- polr(outcome_final ~ Peer   +  PR_char  + Peer*PR_char + 
              startyear, data = mdat, Hess = TRUE)


ggpredict(m_2, term = c("Peer", "PR_char"), 
          condition = c(startyear = mean(mdat$startyear)), 
          vcov.fun = "vcovCL", vcov.type = "HC1", 
          vcov.args = list(cluster = mdat$subclass))

m_2 <- coeftest(m_2, vcov=vcovCL(m_2, factor(mdat$NCP)))
m_2



m_3 <- polr(outcome_final ~ Peer   +  PR_char  + treaty + NGOnumbers + treaty + Bribery +
              startyear + Peer*PR_char, data = mdat, Hess = TRUE)

m_3 <- coeftest(m_3, vcov=vcovCL(m_3, factor(mdat$NCP)))
m_3


m_4 <- polr(outcome_final ~ Peer   +  PR_char + NGOnumbers + treaty + Bribery + legparty + 
              startyear + Peer*PR_char, data = mdat, Hess = TRUE)

m_4 <- coeftest(m_4, vcov=vcovCL(m_4, factor(mdat$NCP)))
m_4

# Table 3 in the Appendix C. Matching Results
stargazer(m_1, m_2, m_3, m_4, digits=2, style = "IO")



############# Domestic NGOs as treatment in matching

df = na.omit(subset(dat, select = c(outcome_final, NGOnumbers, peerrev, pr, incumpower, opppower, treaty, binarypta,
                                    polipush, legparty, NCP_GDPpercap,
                                    provision, provision2, startyear, endyear, Labor, Environment, Bribery, NCP
)))


df$domestic <- NA

for(i in 1:nrow(df)){
  if(df$polipush[i] > 0){df$domestic[i] <- 1}
  else{df$domestic[i] <- 0}
}


m.out1 <- matchit(domestic ~  NCP_GDPpercap + legparty, data = df,
                  method = "cem", cutpoints = list(NCP_GDPpercap = c(39707,50097)), distance = "glm")

summary(m.out1)
mdat2 <- match.data(m.out1)

# Figure 1 in Appendix B. Matching Process
plot(summary(m.out1), xlim=c(0,2))


mdat2$PR_char <- NA

for(i in 1:nrow(mdat2)){
  if(mdat2$pr[i] == 1){mdat2$PR_char[i] <- "PR"}
  else{mdat2$PR_char[i] <- "Non-PR"}
}

mdat2$PR_char <- as.factor(mdat2$PR_char)

mdat2$Domestic <- NA

for(i in 1:nrow(mdat2)){
  if(mdat2$domestic[i] == 1){mdat2$Domestic[i] <- "Home NGOs"}
  else{mdat2$Domestic[i] <- "Foreign"}
}

mdat2$Domestic <- as.factor(mdat2$Domestic)


summary(m.out1, un = FALSE)
plot(summary(m.out1), xlim=c(0,2))

mat_0 <- polr(outcome_final ~  Domestic +  PR_char  + startyear, data = mdat2, Hess = TRUE)
mat_0 <- coeftest(mat_0, vcov=vcovCL(mat_0, factor(mdat2$NCP)))
mat_0


mat_1 <- polr(outcome_final ~  Domestic +  PR_char  + Domestic*PR_char  + startyear, data = mdat2, Hess = TRUE)

ggpredict(mat_1, term = c("Domestic", "PR_char"), 
          condition = c(startyear = mean(mdat$startyear)),
          vcov.fun = "vcovCL", vcov.type = "HC1", 
          vcov.args = list(cluster = mdat2$subclass))

mat_1 <- coeftest(mat_1, vcov=vcovCL(mat_1, factor(mdat2$NCP)))
mat_1



mat_2 <- polr(outcome_final ~  Domestic +  PR_char  + Domestic*PR_char  + treaty + 
                Bribery + startyear, data = mdat2, Hess = TRUE)
mat_2 <- coeftest(mat_2, vcov=vcovCL(mat_2, factor(mdat2$NCP)))
mat_2




mat_3 <- polr(outcome_final ~   Domestic +  PR_char  + log(NCP_GDPpercap)  + treaty + Bribery +
                NGOnumbers + legparty+ startyear + Domestic*PR_char, data = mdat2, Hess = TRUE)

mat_3 <- coeftest(mat_3, vcov=vcovCL(mat_3, factor(mdat2$NCP)))
mat_3

# Table 4 in Appendix C. Matching results
stargazer(mat_0, mat_1, mat_2, mat_3, digits=2, style = "IO")


###########################################################
########################## PR alternative coding with Mixed Systems using Golder & Bormann (Appendix D)
###########################################################

g1 <- polr(outcome_final ~  peerrev + pr3 + startyear, data = dat, Hess = TRUE)
g1 <- coeftest(g1, vcov=vcovCL(g1, factor(dat$NCP)))
g1

g2 <- polr(outcome_final ~  peerrev + pr3  + incumpower + opppower + provision + provision2 + startyear, data = dat, Hess = TRUE)
g2 <- coeftest(g2, vcov=vcovCL(g2, factor(dat$NCP)))
g2


g3 <- polr(outcome_final ~  peerrev + pr3 + NGOnumbers + incumpower + opppower + provision + provision2 + peerrev*pr3 + startyear, data = dat, Hess = TRUE)
g3 <- coeftest(g3, vcov=vcovCL(g3, factor(dat$NCP)))
g3

# Table 5 in Appendix D
stargazer(g1, g2, g3, digits=2, style ="IO")


dat$domestic <- NA

for(i in 1:nrow(dat)){
  if(is.na(dat$polipush[i])==1){dat$domestic[i] <- NA}
  else if(dat$polipush[i] > 0){dat$domestic[i] <- 1}
  else{dat$domestic[i] <- 0}
}


gd1 <- polr(outcome_final ~  domestic + pr3  + startyear, data = dat, Hess = TRUE)
gd1 <- coeftest(gd1, vcov=vcovCL(gd1, factor(dat$NCP)))
gd1

gd2 <- polr(outcome_final ~  domestic + pr3 + incumpower + opppower + provision + provision2 + startyear, data = dat, Hess = TRUE)
gd2 <- coeftest(gd2, vcov=vcovCL(gd2, factor(dat$NCP)))
gd2


gd3 <- polr(outcome_final ~  domestic + pr3 + incumpower + opppower + provision + provision2 + domestic*pr3 + startyear, 
            data = dat, Hess = TRUE)
gd3 <- coeftest(gd3, vcov=vcovCL(gd3, factor(dat$NCP)))
gd3

# Table 6 in Appendix D
stargazer(gd1, gd2, gd3, digits=2, style ="IO")




dat_fig <- aggregate(count ~ pr3 + Decision, data = dat, sum)


# Figure 2 in Appendix D

dat_fig$labely[1] <- 15 + 27 + (55/2)     # Mixed - reject
dat_fig$labely[2]<- 3 + 22 + (29/2)                  # Plural - reject
dat_fig$labely[3] <-  14 + 35 + (43/2)                 # PR - reject
dat_fig$labely[4]<- 15 + (27/2)          # Neutral -Mixed
dat_fig$labely[5]<- 3 + (22/2)          # Neutral -Plural
dat_fig$labely[6]<-   14 + (35/2)    # Neutral -PR
dat_fig$labely[7]<- 15/2         # Recommend - Mixed
dat_fig$labely[8]<- 3/2  #  Recommend - Plurality
dat_fig$labely[9]<- 14/2 # recommend -  PR

ggplot(dat_fig, aes(x = pr3, y = count, fill = Decision)) +
  geom_col(colour = "black") +
  scale_fill_brewer(palette = "Pastel1") +
  geom_text(aes(label = count, y = labely), size = 3) + theme_classic() + 
  xlab("") + ylab("Count of cases")  +  scale_fill_manual(values=c('darkred','lightgray', 'darkgreen'))



#####################################################
################# District size (Appendix J, Figure 4 in the manuscript)
#####################################################

district1 <- polr(outcome_final ~  
                    peerrev + polipush + mdmhnew +  incumpower  + opppower  + startyear, data = dat, Hess = TRUE)

district1  <- coeftest(district1 , vcov=vcovCL(district1 , factor(dat$NCP)))
district1



district2 <- polr(outcome_final ~  
                    peerrev + polipush + mdmhnew +  incumpower  + opppower  + provision + provision2 + NGOnumbers +  startyear, data = dat, Hess = TRUE)

district2  <- coeftest(district2 , vcov=vcovCL(district2 , factor(dat$NCP)))
district2


district3 <- polr(outcome_final ~  
                    peerrev + polipush + mdmhnew +  incumpower  + opppower  + provision + provision2 + treaty + NGOnumbers + 
                    Labor + Environment + Bribery +
                    startyear, data = dat, Hess = TRUE)

# Figure 4 in the manuscript
plot(ggpredict(district3, term = c("mdmhnew"),              
               condition = c(startyear = 2017,
                             provision = mean(dat$provision, na.rm=T), 
                             incumpower = mean(dat$incumpower, na.rm=T),
                             opppower = mean(dat$opppower, na.rm=T), 
                             NGOnumbers = mean(dat$NGOnumbers), 
                             Labor = 0, Environment = 0, Bribery = 0,
                             polipush =1,  peerrev= 1,
                             provision2 = mean(dat$provision, na.rm=T)*mean(dat$provision, na.rm=T)), 
               vcov.fun = "vcovCL", vcov.type = "HC1", vcov.args = list(cluster = dat$subclass))) + 
  xlab("District Magnitude") + 
  ylab("Predicted Probabilities") + 
  ggtitle("") + theme_classic() + scale_color_grey(start=0.7, end = 0.2)


ggpredict(district3, term = c("mdmhnew[1:8]"),              
          condition = c(startyear = 2017,
                        provision = mean(dat$provision, na.rm=T), 
                        incumpower = mean(dat$incumpower, na.rm=T),
                        opppower = mean(dat$opppower, na.rm=T), 
                        NGOnumbers = mean(dat$NGOnumbers), 
                        Labor = 0, Environment = 0, Bribery = 0,
                        polipush =1,  peerrev= 1,
                        provision2 = mean(dat$provision, na.rm=T)*mean(dat$provision, na.rm=T)), 
          vcov.fun = "vcovCL", vcov.type = "HC1", vcov.args = list(cluster = dat$subclass))

district3  <- coeftest(district3 , vcov=vcovCL(district3 , factor(dat$NCP)))
district3

# Table 12 in Appendix J.
stargazer(district1, district2, district3, digits=2, style="IO")



######################################################
############ Cross-validation (Appendix I)
######################################################
# Peer Review * PR
peer_cross <- polr(outcome_final ~  
                     peerrev + pr  + NGOnumbers + incumpower + opppower +  treaty  + binarypta +
                     provision + provision2 + pr*peerrev  + Labor + Environment + Bribery + startyear, data = dat, Hess = TRUE)

peer_cross <- coeftest(peer_cross, vcov=vcovCL(peer_cross, factor(dat$NCP)))
peer_cross


country_list <- unique(dat$NCP)
country_list

cv_holder <- data.frame(country_dropped = country_list,
                        peerrev_pr_int_coef = NA,
                        peerrev_pr_se = NA
)

# Create the for loop and iterate
for (i in 1:nrow(cv_holder)) {
  
  drop_country <- cv_holder$country_dropped[i]
  
  # Drop a country from the dataset
  d_loop <- dat %>% filter(NCP != drop_country)
  
  # Fit the model
  
  fit_loop <- polr(outcome_final ~  
                     peerrev + pr  + NGOnumbers + incumpower + opppower +  treaty  + binarypta +
                     provision + provision2 + pr*peerrev  + Labor + Environment + Bribery + startyear,
                   data = d_loop, # Update the data to be the current subset
                   Hess = TRUE)
  # cluster SEs
  fit_loop <- coeftest(fit_loop, vcov=vcovCL(fit_loop, factor(d_loop$NCP)))
  fit_loop
  
  # Extract and save coefficients / stats from the model
  fit_tidy <- tidy(fit_loop)
  
  
  # peerrev:pr
  cv_holder$peerrev_pr_int_coef[i] <- fit_tidy$estimate[14]
  cv_holder$peerrev_pr_se[i] <- fit_tidy$std.error[14]
  
}

full_sample_row <- data.frame(
  country_dropped = "All Sample",  # You can change the label as needed
  peerrev_pr_int_coef = tidy(peer_cross)$estimate[14],
  peerrev_pr_se = tidy(peer_cross)$std.error[14]
)

cv_holder <- bind_rows(cv_holder, full_sample_row)

# Figure 3 left panel in Appendix I.
cv_holder %>%
  mutate(highlight = ifelse(country_dropped == "All Sample", "Yes", "No")) %>%
  ggplot(
    aes(x = country_dropped, y = peerrev_pr_int_coef)
  ) +
  geom_point() +
  geom_hline(yintercept = tidy(peer_cross)$estimate[14], col = 'red') +
  geom_errorbar(
    aes(ymin = peerrev_pr_int_coef - 1.96 * peerrev_pr_se, ymax = peerrev_pr_int_coef + 1.96 * peerrev_pr_se),
    width = 0.2) +
  scale_color_manual(values = c("No" = "black", "Yes" = "red")) +  # Set colors for points
  theme_clean() +
  labs(
    x = "",
    y = "Coefficient Estimate on Interaction Term",
    title = "Coefficient Estimates on Peer Review:PR After Dropping Each Country",
    subtitle = "Coefficient estimate from the full sample plotted in red"
  ) +
  scale_x_discrete(guide = guide_axis(angle = 90)) + ylim(0.8,1.8)


# Domestic NGO * PR
domestic_cross <- polr(outcome_final ~  
                         polipush + pr  + NGOnumbers + incumpower + opppower +  treaty  + binarypta +
                         provision + provision2 + pr*polipush  + peerrev + Labor + Environment + Bribery + startyear, data = dat, Hess = TRUE)

domestic_cross <- coeftest(domestic_cross, vcov=vcovCL(domestic_cross, factor(dat$NCP)))
domestic_cross


country_list <- unique(dat$NCP)
country_list

cv_holder <- data.frame(country_dropped = country_list,
                        domestic_pr_int_coef = NA,
                        domestic_pr_se = NA
)

# Create the for loop and iterate
for (i in 1:nrow(cv_holder)) {
  
  drop_country <- cv_holder$country_dropped[i]
  
  # Drop a country from the dataset
  d_loop <- dat %>% filter(NCP != drop_country)
  
  # Fit the model
  
  fit_loop <- polr(outcome_final ~  
                     polipush + pr  + NGOnumbers + incumpower + opppower +  treaty  + binarypta +
                     provision + provision2 + pr*polipush  + peerrev + Labor + 
                     Environment + Bribery + startyear, data = d_loop, Hess = TRUE)
  # cluster SEs
  fit_loop <- coeftest(fit_loop, vcov=vcovCL(fit_loop, factor(d_loop$NCP)))
  fit_loop
  
  # Extract and save coefficients / stats from the model
  fit_tidy <- tidy(fit_loop)
  
  # domestic:pr
  cv_holder$domestic_pr_int_coef[i] <- fit_tidy$estimate[15]
  cv_holder$domestic_pr_se[i] <- fit_tidy$std.error[15]
  
}

full_sample_row <- data.frame(
  country_dropped = "All Sample",  
  domestic_pr_int_coef = tidy(domestic_cross)$estimate[15],
  domestic_pr_se = tidy(domestic_cross)$std.error[15]
)

cv_holder <- bind_rows(cv_holder, full_sample_row)


# Figure 3 right panel in Appendix I.
cv_holder %>%
  mutate(highlight = ifelse(country_dropped == "All Sample", "Yes", "No")) %>%
  ggplot(
    aes(x = country_dropped, y = domestic_pr_int_coef)
  ) +
  geom_point() +
  geom_hline(yintercept = tidy(domestic_cross)$estimate[15], col = 'red') +
  geom_errorbar(
    aes(ymin = domestic_pr_int_coef - 1.96 * domestic_pr_se, ymax = domestic_pr_int_coef + 1.96 * domestic_pr_se),
    width = 0.2) +
  scale_color_manual(values = c("No" = "black", "Yes" = "red")) +  # Set colors for points
  theme_clean() +
  labs(
    x = "",
    y = "Coefficient Estimate on Interaction Term",
    title = "Coefficient Estimates on domestic NGO:PR After Dropping Each Country",
    subtitle = "Coefficient estimate from the full sample plotted in red"
  ) +
  scale_x_discrete(guide = guide_axis(angle = 90)) 



######################################################
########################## Figure 2 in the manuscript
######################################################

d_plot <- dat %>%
  group_by(NCP, outcome_order) %>%
  summarise(Outcome = n()) %>%  ungroup()

d_plot %>% print(n=100)


d_plot$outcome_order <- as.numeric(d_plot$outcome_order)

for(i in 1:nrow(d_plot)){
  if(d_plot$outcome_order[i] == 1){d_plot$outcome_order[i] <- "1.Vindication"}
  else if(d_plot$outcome_order[i] == 2){d_plot$outcome_order[i] <- "2.Neutral"}
  else if(d_plot$outcome_order[i] == 3){d_plot$outcome_order[i] <- "3.Recommendation"}
}


d_plot$Outcome <- as.numeric(d_plot$Outcome)
d_plot$NCP <- as.character(d_plot$NCP)

d_plot %>% print(n = 100)

sum(dat$outcome_reject)
sum(dat[dat$outcome_order == 3,]$count)/243

d_plot$NCP <- factor(d_plot$NCP, 
                     levels = c("Poland", "Peru", "Austria", "Chile",
                                "Czech Republic", "Israel", "Luxembourg", "Ireland",
                                "Mexico", "Finland", "Spain", "Australia",
                                "Sweden","Brazil", "Italy",  "France", "New Zealand", "Japan", "Korea",
                                "Switzerland", "Denmark", "Norway", "Belgium", "Canada",
                                "Netherlands", "USA", "Germany", "UK"))


ggplot(d_plot, aes(x = NCP, y = Outcome, fill = outcome_order)) +
  geom_bar(stat="identity") + coord_flip() + theme_bw() +
  scale_fill_manual(values=c("black", "grey", "red4"), name = "Decision") + ylab("Count of Cases") + xlab("Government NCP")+
  theme(legend.position = "bottom", legend.direction = "horizontal") 



######################################################
########################## EU as a variable (Appendix L)
######################################################

eu1 <- polr(outcome_final ~  
              peerrev  + pr + NGOnumbers + incumpower  + opppower +
              provision + provision + provision2 + eucountry + 
              startyear, data = dat, Hess = TRUE)
eu1   <- coeftest(eu1  , vcov=vcovCL(eu1  , factor(dat$NCP)))
eu1


eu2 <- polr(outcome_final ~  
              peerrev  + pr + NGOnumbers + incumpower  + opppower + treaty + 
              provision + provision + provision2 + eucountry + Labor + Environment + Bribery +
              startyear, data = dat, Hess = TRUE)
eu2   <- coeftest(eu2  , vcov=vcovCL(eu2  , factor(dat$NCP)))
eu2

eu3 <- polr(outcome_final ~  
              peerrev  + pr + NGOnumbers + incumpower  + opppower + treaty + 
              provision + provision + provision2 + eucountry + Labor + Environment + Bribery + peerrev*pr +
              startyear, data = dat, Hess = TRUE)
eu3   <- coeftest(eu3  , vcov=vcovCL(eu3  , factor(dat$NCP)))
eu3




eu1_ngo <- polr(outcome_final ~  
                  domestic  + pr + NGOnumbers + incumpower  + opppower +
                  provision + provision + provision2 + eucountry + 
                  startyear, data = dat, Hess = TRUE)
eu1_ngo   <- coeftest(eu1_ngo  , vcov=vcovCL(eu1_ngo  , factor(dat$NCP)))
eu1_ngo


eu2_ngo <- polr(outcome_final ~  
                  domestic  + pr + NGOnumbers + incumpower  + opppower + treaty + 
                  provision + provision + provision2 + eucountry + Labor + Environment + Bribery +
                  startyear, data = dat, Hess = TRUE)
eu2_ngo   <- coeftest(eu2_ngo  , vcov=vcovCL(eu2_ngo  , factor(dat$NCP)))
eu2_ngo


eu3_ngo <- polr(outcome_final ~  
                  domestic  + pr + NGOnumbers + incumpower  + opppower + treaty + 
                  provision + provision + provision2 + eucountry + Labor + Environment + Bribery + domestic*pr +
                  startyear, data = dat, Hess = TRUE)
eu3_ngo   <- coeftest(eu3_ngo  , vcov=vcovCL(eu3_ngo  , factor(dat$NCP)))
eu3_ngo

# Table 14 in Appendix L
stargazer(eu1, eu2, eu3, eu1_ngo, eu2_ngo, eu3_ngo, digits=2, style = "IO")


#######################################################
########## Green party (Appendix M)
#######################################################

eco1 <- polr(outcome_final ~ peerrev + domestic + ecoseat + incumpower + opppower + 
               startyear, data = dat, Hess = TRUE)

ggpredict(eco1, term = c("ecoseat"), 
          condition = c(startyear = mean(dat$startyear),
                        provision = mean(dat$provision, na.rm=T), incumpower = mean(dat$incumpower, na.rm=T), 
                        opppower = mean(dat$opppower, na.rm=T),
                        peerrev = 0, domestic = 1), 
          vcov.fun = "vcovCL", vcov.type = "HC1", 
          vcov.args = list(cluster = dat$subclass))


eco1 <- coeftest(eco1, vcov=vcovCL(eco1, factor(dat$NCP)))
eco1



eco2 <- polr(outcome_final ~ peerrev + domestic + ecoseat + incumpower + opppower + NGOnumbers +
               provision + provision2 + 
               startyear, data = dat, Hess = TRUE)

# this is Figure 8 in Appendix 
plot(ggpredict(eco2, term = c("ecoseat"), 
               condition = c(startyear = mean(dat$startyear),
                             provision = mean(dat$provision, na.rm=T), incumpower = mean(dat$incumpower, na.rm=T), 
                             opppower = mean(dat$opppower, na.rm=T),
                             NGOnumbers = mean(dat$NGOnumbers),
                             provision2 = mean(dat$provision, na.rm=T)*mean(dat$provision, na.rm=T)), 
               vcov.fun = "vcovCL", vcov.type = "HC1", 
               vcov.args = list(cluster = dat$subclass)), connect.lines = TRUE, ci.style = "dot") + xlab("") + theme_classic() + 
  ggtitle("") + xlab("Green Party Seat Share") + ylab("Predicted Probabilites of \n Vindication, Neutral, Recommendation")


eco2 <- coeftest(eco2, vcov=vcovCL(eco2, factor(dat$NCP)))
eco2


eco3 <- polr(outcome_final ~ peerrev + domestic + ecoseat + incumpower + opppower + NGOnumbers +
               provision + provision2 + treaty + Labor + Environment + Bribery + 
               startyear, data = dat, Hess = TRUE)

ggpredict(eco3, term = c("ecoseat[0,0.13]"), 
          condition = c(startyear = mean(dat$startyear),
                        provision = mean(dat$provision, na.rm=T), incumpower = mean(dat$incumpower, na.rm=T), 
                        opppower = mean(dat$opppower, na.rm=T),
                        NGOnumbers = mean(dat$NGOnumbers),
                        provision2 = mean(dat$provision, na.rm=T)*mean(dat$provision, na.rm=T)), 
          vcov.fun = "vcovCL", vcov.type = "HC1", 
          vcov.args = list(cluster = dat$subclass))


eco3 <- coeftest(eco3, vcov=vcovCL(eco3, factor(dat$NCP)))
eco3

# Table 15 in Appendix M
stargazer(eco1, eco2, eco3, digits=2, style="IO")


####################################
################## Neoliberalism & Corporatism (Appendix H)
####################################

neo1 <- polr(outcome_final ~  overall + 
               peerrev + pr + incumpower + opppower + 
               provision + provision2  + peerrev*pr + NGOnumbers + 
               Labor + Environment + Bribery + startyear, data = dat, Hess = TRUE)

neo1  <- coeftest(neo1 , vcov=vcovCL(neo1 , factor(dat$NCP)))
neo1 

neo2 <- polr(outcome_final ~  tradefreedom + 
               peerrev + pr + incumpower + opppower + 
               provision + provision2  + peerrev*pr + NGOnumbers + 
               Labor + Environment + Bribery + startyear, data = dat, Hess = TRUE)

neo2  <- coeftest(neo2 , vcov=vcovCL(neo2 , factor(dat$NCP)))
neo2 


neo3 <- polr(outcome_final ~  overall + 
               polipush + pr + incumpower + opppower + 
               provision + provision2  + polipush*pr + NGOnumbers + peerrev +
               Labor + Environment + Bribery + startyear, data = dat, Hess = TRUE)

neo3  <- coeftest(neo3 , vcov=vcovCL(neo3, factor(dat$NCP)))
neo3

neo4 <- polr(outcome_final ~  tradefreedom + 
               polipush + pr + incumpower + opppower + 
               provision + provision2  + polipush*pr + NGOnumbers + peerrev + 
               Labor + Environment + Bribery + startyear, data = dat, Hess = TRUE)
neo4  <- coeftest(neo4 , vcov=vcovCL(neo4 , factor(dat$NCP)))
neo4


cor1 <- polr(outcome_final ~  CorpAll + 
               peerrev + pr + incumpower + opppower + 
               provision + provision2  + peerrev*pr + NGOnumbers + 
               Labor + Environment + Bribery + startyear, data = dat, Hess = TRUE)

cor1  <- coeftest(cor1 , vcov=vcovCL(cor1 , factor(dat$NCP)))
cor1

cor2 <- polr(outcome_final ~  CorpAll + 
               peerrev + polipush + pr + incumpower + opppower + 
               provision + provision2  + polipush*pr + NGOnumbers + 
               Labor + Environment + Bribery + startyear, data = dat, Hess = TRUE)

cor2  <- coeftest(cor2 , vcov=vcovCL(cor2 , factor(dat$NCP)))
cor2

# Table 11 in Appendix H
stargazer(neo1, neo2, neo3, neo4, cor1, cor2, digits=2, style = "IO")


####################################
################## Campaign finance (Appendix N)
####################################


hummel1 <- polr(outcome_final ~ domestic + peerrev + ecoseat + incumpower + opppower +  PFSI   + Labor + Environment + Bribery + provision + provision2 +  
                  startyear, data = dat, Hess = TRUE)

hummel1 <- coeftest(hummel1, vcov=vcovCL(hummel1, factor(dat$NCP)))
hummel1


hummel2 <- polr(outcome_final ~ domestic + peerrev + mdmhnew + incumpower + opppower +  PFSI   + Labor + Environment + Bribery + provision + provision2 + 
                  startyear, data = dat, Hess = TRUE)
hummel2 <- coeftest(hummel2, vcov=vcovCL(hummel2, factor(dat$NCP)))
hummel2



vdem1 <- polr(outcome_final ~ domestic + peerrev + mdmhnew + incumpower + opppower +  pargroup   + Labor + Environment + Bribery + provision + provision2 + 
                startyear, data = dat, Hess = TRUE)
vdem1 <- coeftest(vdem1, vcov=vcovCL(vdem1, factor(dat$NCP)))
vdem1


vdem2 <- polr(outcome_final ~ domestic + peerrev + ecoseat + incumpower + opppower +  pargroup   + Labor + Environment + Bribery  + provision + provision2 + 
                startyear, data = dat, Hess = TRUE)
vdem2 <- coeftest(vdem2, vcov=vcovCL(vdem2, factor(dat$NCP)))
vdem2

vdem3 <- polr(outcome_final ~ peerrev + PR +  PR*peerrev + incumpower + opppower +  pargroup   + Labor + Environment + Bribery  + provision + provision2 + 
                startyear, data = dat, Hess = TRUE)

vdem3 <- coeftest(vdem3, vcov=vcovCL(vdem3, factor(dat$NCP)))
vdem3


vdem4 <- polr(outcome_final ~ domestic + PR +  domestic*PR + incumpower + opppower +  pargroup   + Labor + Environment + Bribery   + provision + provision2 + 
                startyear, data = dat, Hess = TRUE)

vdem4 <- coeftest(vdem4, vcov=vcovCL(vdem4, factor(dat$NCP)))
vdem4

vdem5 <- polr(outcome_final ~ domestic + peerrev + mdmhnew + incumpower + opppower +  parfund   + Labor + Environment + Bribery + provision + provision2 + 
                startyear, data = dat, Hess = TRUE)
vdem5 <- coeftest(vdem5, vcov=vcovCL(vdem5, factor(dat$NCP)))
vdem5


vdem6 <- polr(outcome_final ~ domestic + peerrev + ecoseat + incumpower + opppower +  parfund   + Labor + Environment + Bribery  + provision + provision2 + 
                startyear, data = dat, Hess = TRUE)
vdem6 <- coeftest(vdem6, vcov=vcovCL(vdem6, factor(dat$NCP)))
vdem6

vdem7 <- polr(outcome_final ~ peerrev + PR +  PR*peerrev + incumpower + opppower +  parfund   + Labor + Environment + Bribery  + provision + provision2 + 
                startyear, data = dat, Hess = TRUE)
vdem7 <- coeftest(vdem7, vcov=vcovCL(vdem7, factor(dat$NCP)))
vdem7


vdem8 <- polr(outcome_final ~ domestic + PR +  domestic*PR + incumpower + opppower +  parfund   + Labor + Environment + Bribery   + provision + provision2 + 
                startyear, data = dat, Hess = TRUE)
vdem8 <- coeftest(vdem8, vcov=vcovCL(vdem8, factor(dat$NCP)))
vdem8

# Table 16 in Appendix N
stargazer(hummel1, hummel2, vdem2,vdem1, vdem3, vdem4,  digits=2, style = "IO")



####################################
#################### time trend (Not in the Supplementary Information or manuscript)
####################################

peer_cubic <- polr(outcome_final ~  
                     peerrev + pr  + NGOnumbers + incumpower + opppower +  treaty  + binarypta +
                     provision + provision2 + pr*peerrev  + Labor + Environment + Bribery + 
                     startyearnew + startyearnew2 + startyearnew3, data = dat, Hess = TRUE)

peer_cubic <- coeftest(peer_cubic, vcov=vcovCL(peer_cubic, factor(dat$NCP)))
peer_cubic


domestic_cubic <- polr(outcome_final ~  
                         polipush + pr + NGOnumbers + incumpower + opppower +
                         provision + provision2 +  polipush*pr +  treaty + binarypta + 
                         Labor + Environment + Bribery + 
                         startyearnew + startyearnew2 + startyearnew3 , data = dat, Hess = TRUE)
domestic_cubic <- coeftest(domestic_cubic, vcov=vcovCL(domestic_cubic, factor(dat$NCP)))
domestic_cubic


eco_cubic <- polr(outcome_final ~ peerrev  + domestic + ecoseat + incumpower + opppower + provision2 + provision + 
                    startyearnew + startyearnew2 + startyearnew3, data = dat, Hess = TRUE)
eco_cubic <- coeftest(eco_cubic, vcov=vcovCL(eco_cubic, factor(dat$NCP)))
eco_cubic


district_cubic <- polr(outcome_final ~  
                         peerrev + domestic + mdmhnew +  incumpower  + opppower  + 
                         provision + provision2 + treaty + NGOnumbers + 
                         Labor + Environment + Bribery +
                         startyearnew + startyearnew2 + startyearnew3, data = dat, Hess = TRUE)

district_cubic <- coeftest(district_cubic, vcov=vcovCL(district_cubic, factor(dat$NCP)))
district_cubic



peer_tfe <- polr(outcome_final ~  
                   peerrev + pr   + incumpower + opppower +
                   provision + provision2 + peerrev*pr   +
                   factor(startyear)  , data = dat, Hess = TRUE)

peer_tfe <- coeftest(peer_tfe, vcov=vcovCL(peer_tfe, factor(dat$NCP)))
peer_tfe


domestic_tfe <- polr(outcome_final ~  
                       polipush + pr + incumpower + opppower  +
                       provision + provision2 +  polipush*pr  +
                       factor(startyear) , data = dat, Hess = TRUE)

domestic_tfe <- coeftest(domestic_tfe, vcov=vcovCL(domestic_tfe, factor(dat$NCP)))
domestic_tfe


eco_tfe <- polr(outcome_final ~ peerrev  + domestic + ecoseat + incumpower + opppower + 
                  factor(startyear), data = dat, Hess = TRUE)
eco_tfe <- coeftest(eco_tfe, vcov=vcovCL(eco_tfe, factor(dat$NCP)))
eco_tfe


district_tfe <- polr(outcome_final ~  
                       peerrev + domestic + mdmhnew +  incumpower  + opppower  + 
                       provision + provision2 + 
                       factor(startyear), data = dat, Hess = TRUE)

district_tfe <- coeftest(district_tfe, vcov=vcovCL(district_tfe, factor(dat$NCP)))
district_tfe


######################################################
############ alternative coding of the DV (Appendix P)
######################################################


dat$outcome_four <- as.factor(dat$outcome_four)

alt1 <- polr(outcome_four ~  
               peerrev + pr + incumpower  + opppower  + treaty + peerrev + polipush +
               provision + provision2 + pr*peerrev + Labor + Environment + Bribery +
               startyear, data = dat, Hess = TRUE)

alt1 <- coeftest(alt1, vcov=vcovCL(alt1, factor(dat$NCP)))
alt1



alt2 <- polr(outcome_four ~  
               polipush + pr + incumpower  + opppower  + treaty + peerrev +
               provision + provision2 + pr*polipush + Labor + Environment + Bribery +
               startyear, data = dat, Hess = TRUE)

alt2 <- coeftest(alt2, vcov=vcovCL(alt2, factor(dat$NCP)))
alt2


alt3 <- polr(outcome_four ~  
               polipush +  + incumpower  + opppower  + treaty + peerrev + ecoseat +
               provision + provision2  + Labor + Environment + Bribery +
               startyear, data = dat, Hess = TRUE)

alt3 <- coeftest(alt3, vcov=vcovCL(alt3, factor(dat$NCP)))
alt3


alt4 <- polr(outcome_four ~  
               polipush +  + incumpower  + opppower  + treaty + peerrev + mdmhnew +
               provision + provision2  + Labor + Environment + Bribery +
               startyear, data = dat, Hess = TRUE)

alt4 <- coeftest(alt4, vcov=vcovCL(alt4, factor(dat$NCP)))
alt4

# Table 18 in Appenedix P
stargazer(alt1, alt2, alt3, alt4, digits=2, style= "IO")

###################################### 
################### crosstabs in Appendices A and Q
###################################### 

# Table 2 in Appendix A
dat %>%
  group_by(outcome_final, pr) %>%
  tally() %>%
  spread(outcome_final, n)

# Table 19 in Appendix Q
dat %>%
  group_by(pr, peerrev) %>%
  tally() %>%
  spread(pr, n)

# Table 20 in Appendix Q
dat %>%
  group_by(pr, peerrev_1yr) %>%
  tally() %>%
  spread(pr, n)


m_pralt <- polr(outcome_final ~  
                  peerrev_1yr + pr   + NGOnumbers + incumpower + opppower +  treaty + binarypta + 
                  provision + provision2 + pr*peerrev_1yr  + Labor + Environment + Bribery + 
                  startyear, data = dat, Hess = TRUE)

m_pralt <- coeftest(m_pralt, vcov=vcovCL(m_pralt, factor(dat$NCP)))
m_pralt

# Table 21 in Appendix Q
stargazer(m_pralt, digits=2, style = "IO")



########################################################
######## PTA Mechanism Appendix F
#######################################################

# peer review * pr 
m_4 <- polr(outcome_final ~  
              peerrev + pr  + NGOnumbers + incumpower + opppower +  treaty  + binarypta +
              provision + provision2 + pr*peerrev  + Labor + Environment + Bribery + startyear, data = dat, Hess = TRUE)

m_4 <- coeftest(m_4, vcov=vcovCL(m_4, factor(dat$NCP)))
m_4

# ESR 
esr_pr <- polr(outcome_final ~  
                 peerrev + pr  + NGOnumbers + incumpower + opppower +  treaty  + esr_all_lta +
                 provision + provision2 + pr*peerrev  + Labor + Environment + Bribery + startyear, data = dat, Hess = TRUE)

esr_pr <- coeftest(esr_pr, vcov=vcovCL(esr_pr, factor(dat$NCP)))
esr_pr


# EP 
ep_pr <- polr(outcome_final ~  
                peerrev + pr  + NGOnumbers + incumpower + opppower +  treaty  + ep_all_lta +
                provision + provision2 + pr*peerrev  + Labor + Environment + Bribery + startyear, data = dat, Hess = TRUE)

ep_pr <- coeftest(ep_pr, vcov=vcovCL(ep_pr, factor(dat$NCP)))
ep_pr



# CP
cp_pr <- polr(outcome_final ~  
                peerrev + pr  + NGOnumbers + incumpower + opppower +  treaty  + cpr_all_lta +
                provision + provision2 + pr*peerrev  + Labor + Environment + Bribery + startyear, data = dat, Hess = TRUE)

cp_pr <- coeftest(cp_pr, vcov=vcovCL(cp_pr, factor(dat$NCP)))
cp_pr

# Table 9 in Appendix F
stargazer(m_4, esr_pr, ep_pr, cp_pr, digits=2, style = "IO")


# domestic * pr 

d4 <- polr(outcome_final ~  
             polipush + pr + NGOnumbers + incumpower + opppower +
             provision + provision2 +  polipush*pr +  treaty + binarypta + peerrev +
             Labor + Environment + Bribery + startyear , data = dat, Hess = TRUE)

summary(d4)
d4 <- coeftest(d4, vcov=vcovCL(d4, factor(dat$NCP)))
d4

# esr
esr_d <- polr(outcome_final ~  polipush + pr + NGOnumbers + incumpower + opppower +
                provision + provision2 +  polipush*pr  + esr_all_lta  + peerrev + treaty +
                Labor + Environment + Bribery + startyear  , data = dat, Hess = TRUE)

esr_d <- coeftest(esr_d, vcov=vcovCL(esr_d, factor(dat$NCP)))
esr_d



# ep
ep_d <- polr(outcome_final ~  polipush + pr + NGOnumbers + incumpower + opppower +
               provision + provision2 +  polipush*pr  + ep_all_lta  + peerrev + treaty +
               Labor + Environment + Bribery + startyear  , data = dat, Hess = TRUE)

ep_d <- coeftest(ep_d, vcov=vcovCL(ep_d, factor(dat$NCP)))
ep_d



# cpr
cpr_d <- polr(outcome_final ~  polipush + pr + NGOnumbers + incumpower + opppower +
                provision + provision2 +  polipush*pr  + cpr_all_lta  + peerrev + treaty +
                Labor + Environment + Bribery + startyear  , data = dat, Hess = TRUE)

cpr_d <- coeftest(cpr_d, vcov=vcovCL(cpr_d, factor(dat$NCP)))
cpr_d

# Table 10 in Appendix F
stargazer(d4, esr_d, ep_d, cpr_d, digits=2, style = "IO")


########################################
############ Divided Government (Appendix K)
########################################


divided_house_domestic <- polr(outcome_final ~  polipush + pr + NGOnumbers + incumpower + opppower + oppmajh + 
                                 provision + provision2 +  polipush*pr  + binarypta  + peerrev + treaty +
                                 Labor + Environment + Bribery + startyear  , data = dat, Hess = TRUE)

divided_house_domestic <- coeftest(divided_house_domestic, vcov=vcovCL(divided_house_domestic, factor(dat$NCP)))
divided_house_domestic


divided_house_pr <- polr(outcome_final ~  peerrev + pr + NGOnumbers + incumpower + opppower + oppmajh + polipush +
                           provision + provision2 +  peerrev*pr  + binarypta   + treaty +
                           Labor + Environment + Bribery + startyear  , data = dat, Hess = TRUE)

divided_house_pr <- coeftest(divided_house_pr, vcov=vcovCL(divided_house_pr, factor(dat$NCP)))
divided_house_pr

# Table 13 in Appendix K
stargazer(divided_house_pr, divided_house_domestic, digits=2, style = "IO")



##################################
################# Figure 1 in the manuscript
##################################

setwd("/Users/LEEB15/Library/CloudStorage/Dropbox/OECD_Guideline/data/")
figdat1 <- read.csv("manifesto_uk_plot.csv")

figdat1 <- figdat1 %>%
  filter(Party %in% c("Labour", "Conservative", "Green"))

a1 <- ggplot(figdat1, aes(x=Year, y=Salience, group=Party)) +
  geom_line(aes(colour = Party, linetype = Party))+
  geom_point(aes(shape=Party)) + theme_classic() +scale_color_manual(values=c("blue","green","red")) + 
  scale_linetype_manual(values=c("dotted", "solid", "twodash")) + ggtitle("2. UK") 
a1

figdat2 <- read.csv("manifesto_dutch_plot.csv")
a2 <- ggplot(figdat2, aes(x=Year, y=Salience, group=Party)) +
  geom_line(aes(colour = Party, linetype = Party))+
  geom_point(aes(shape=Party)) + theme_classic() +scale_color_manual(values=c("green","red", "blue")) + 
  scale_linetype_manual(values=c("solid", "twodash", "dotted"))  + ggtitle("3. Netherlands")
a2

figdat3 <- read.csv("manifesto_usa_plot.csv")
a3 <- ggplot(figdat3, aes(x=Year, y=Salience, group=Party)) +
  geom_line(aes(colour = Party, linetype = Party)) +
  geom_point(aes(shape=Party)) + theme_classic() +scale_color_manual(values=c("red", "blue")) + 
  scale_linetype_manual(values=c("twodash", "dotted"))  + ggtitle("1. USA") + ylim(0,3)
a3
